﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "fpts")]
    public static unsafe double fpts(double a, double b, double eps, IntPtr f_x_ptr)
    {
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return fpts(a, b, eps);
    }

    /// <summary>
    /// 自适应梯形求积法
    /// f计算被积函数值f(x)的函数名。
    /// </summary>
    /// <param name="a">积分下限。</param>
    /// <param name="b">积分上限。要求b>a。</param>
    /// <param name="eps">积分精度要求。</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double fpts(double a, double b, double eps)
    {
        double h, t, f0, f1, t0;

        h = b - a;
        t = 0.0;
        f0 = f_x(a);
        f1 = f_x(b);
        t0 = h * (f0 + f1) / 2.0;
        fpts_ppp(a, b, h, f0, f1, t0, eps, &t);
        return (t);
    }

    /// <summary>
    /// 自适应梯形求积法
    /// 递归函数
    /// </summary>
    /// <param name="x0"></param>
    /// <param name="x1"></param>
    /// <param name="h"></param>
    /// <param name="f0"></param>
    /// <param name="f1"></param>
    /// <param name="t0"></param>
    /// <param name="eps"></param>
    /// <param name="t"></param>
    private static unsafe void fpts_ppp(double x0, double x1, double h, double f0, double f1, double t0, double eps, double* t)
    {
        double x, f, t1, t2, p, g, eps1;

        x = x0 + h / 2.0;
        f = f_x(x);
        t1 = h * (f0 + f) / 4.0;
        t2 = h * (f + f1) / 4.0;
        p = Math.Abs(t0 - (t1 + t2));
        if ((p < eps) || (h < eps))
        {
            *t = *t + (t1 + t2);
            return;
        }
        else
        {
            g = h / 2.0;
            eps1 = eps / 1.4;
            fpts_ppp(x0, x, g, f0, f, t1, eps1, t);
            fpts_ppp(x, x1, g, f, f1, t2, eps1, t);
            return;
        }
    }

    /*
    // 自适应梯形求积法例
      int main()
      { 
          double a,b,eps,t,fptsf(double);
          a=-1.0; b=1.0; eps=0.000001; 
          t=fpts(a,b,eps,fptsf);
          cout <<"t = " <<t <<endl;
          return 0;
      }
    // 计算被积函数值
      double fptsf(double x)
      { 
          double y;
          y=1.0/(1.0+25.0*x*x);
          return(y);
      }
    */
}

